##线性时间序列基本概念 
setwd("D:\\Project\\P_Rproject\\data\\time\\ch 2-3")
mibm2697=read.table(file="m-ibm2697.txt")
mvw2697=read.table(file="m-vw2697.txt")
head(mvw2697)
par(mfrow=c(2,1))

acf(mibm2697,main="ACF of IBM simple return")  ##ACF图 
acf(log(1+mibm2697),main="ACF of IBM log return") 

acf(mvw2697,main="ACF of Value-weighted index simple return")
acf(log(1+mvw2697),main="ACF of Value-weighted index log return")


Box.test(mibm2697,lag=5,type="Ljung")  ##simple return  BOX检验 
Box.test(mibm2697,lag=10,type="Ljung")
Box.test(log(1+mibm2697),lag=5,type="Ljung") 
Box.test(log(1+mibm2697),lag=10,type="Ljung") 

Box.test(mvw2697,lag=5,type="Ljung")
Box.test(mvw2697,lag=10,type="Ljung")
Box.test(log(1+mvw2697),lag=5,type="Ljung")
Box.test(log(1+mvw2697),lag=10,type="Ljung")

##############自回归模型 AR（P)
###按照BOX-Jenkins  
##识别 
pvw=pacf(mvw2697,lag=10)  ###求V-W指数的PACF 
pvw
round(pvw$acf,2)   ##保留两位小数

arvw=ar(mvw2697,order=10)
arvw$aic     ###注意此处AIC是经过调整了的 


###V-W指数的AR(3)模型参数估计
mvw2697=ts(mvw2697) ##注意ts函数,转换为时间序列格式
head(mvw2697)
armvw3=ar.ols(mvw2697,order=3,demean=F,intercept=T)
?ar.ols
armvw3
armvw3$asy.se.coef #标准误差
summary(armvw3)
###检验
Box.test(armvw3$resid,lag=10,type="Ljung")
##注意：此处Q统计量服从自由度为（m-p)=7卡方分布
##默认给出的P值不对
1-pchisq(15.8,7)  ###求自由度为7的P值

armvw5=ar.ols(mvw2697,aic=F,order=5,demean=F,intercept=T)
armvw5
Box.test(armvw5$resid,lag=10,type="Ljung")
##注意：此处Q统计量服从自由度为（m-p)=5卡方分布
##默认给出的P值不对
1-pchisq(11.2,5)  ###求自由度为5的P值 
###预测 
##AR(3)预测 
length(mvw2697)
armvw3to858=ar.ols(mvw2697[1:858],aic=F,order=3,demean=F,intercept=T)
vwp3=predict(armvw3to858,n.ahead=6)
vwp3
##AR(5)预测 
armvw5to858=ar.ols(mvw2697[1:858],aic=F,order=5,demean=F,intercept=T)
vwp5=predict(armvw5to858,n.ahead=6) ##向前6步预测
vwp5
mvw2697[859:864] ###真实值 
##画预测效果图 

U=vwp5$pred+1.916*vwp5$se
L=vwp5$pred-1.96*vwp5$se
par(mfrow=c(1,1))
ts.plot(mvw2697,vwp5$pred,col=1:2,xlim=c(840,864))
lines(U,col="blue",lty="dashed")
lines(L,col="blue",lty="dashed")
###############滑动平均模型MA(q）
###MA(q)识别
acf(mvw2697) 
###MA(q)估计与诊断 
fit1=arima(mvw2697,order=c(0,0,9),method="ML") #ma(9)
Box.test(fit1$resid,10)
1-pchisq(0.1552,1)

library(FinTS)
fit2=ARIMA(mvw2697,order=c(0,0,9),method="ML") 
fit2$Box.test
tsdiag(fit2)
###MA(q)预测
fit1=arima(mvw2697[1:854],order=c(0,0,9),method="ML") 
mapre=predict(fit1,n.ahead=3)

##############自回归滑动平均ARMA(p,q),ARIMA(p,d,q) 
###ARMA模拟分析 
ts.sim=arima.sim(list(order = c(1,0,1), ar = 0.9,ma=0.8), n = 200)
par(mfrow=c(3,1))
ts.plot(ts.sim)
acf(ts.sim)
pacf(ts.sim)
#####3M公司ARMA分析 
dat=read.table(file="m-3m4697.txt")
head(dat)
m3m4697=dat[,2]
dates=as.Date(as.character(dat[,1]),format="%Y%m%d") 
plot(dates,m3m4697,type="l")
acf(m3m4697)
pacf(m3m4697)
##利用AIC删选模型,AIC越小越好

arima(m3m4697,order=c(0,0,0),method="ML")
arima(m3m4697,order=c(1,0,0),method="ML")
arima(m3m4697,order=c(2,0,0),method="ML")
arima(m3m4697,order=c(0,0,1),method="ML")
arima(m3m4697,order=c(0,0,2),method="ML")
arima(m3m4697,order=c(1,0,1),method="ML")
arima(m3m4697,order=c(1,0,2),method="ML")
###残差平稳性检验 
arma11=arima(m3m4697,order=c(1,0,1),method="ML")
Box.test(arma11$resid,10,type="Ljung")
1-pchisq(8.0674,8)
####预测
predict(arma11,n.ahead=3)

#################单位根非平稳性 

arima(log(1+m3m4697),order=c(1,0,1),method="ML")
arima(log(1+m3m4697),order=c(0,0,0),method="ML")
mean(log(1+m3m4697))
p=cumsum(log(1+m3m4697))    ##求价格 
pstar=cumsum(log(1+m3m4697)-0.0115)
par(mfrow=c(1,1))
plot(p)
lines(pstar,col="blue")
abline(a=0,b=0.0115,col="red")

arima(p,order=c(0,1,1),method="ML")

###Example 2.2 
dat=read.table(file="q-gdp4703.txt") ####注意工作目录 
gdp=dat[,3]
par(mfrow=c(2,2))
plot(gdp,type="l")
dgdp=diff(gdp)
plot(dgdp,type="l")
acf(gdp)
pacf(gdp)

library(urca)
summary(ur.df(gdp,type="trend",selectlags="AIC"))
summary(ur.kpss(gdp)) ###注意这个原假设是平稳的
###单位根检验更好的函数 
library(fUnitRoots)
help(package="fUnitRoots")
urdfTest(gdp,type="ct")


#####Example 2.2 S&P500
dat=read.table(file="d-sp9003lev.txt")
sp=dat[,2]
plot(sp,type="l")
urdfTest(sp,type="ct")
urdfTest(sp,type="c")
urdfTest(sp,type="nc")



################################季节模型 
jnj=read.table(file="q-jnj.txt")
jnj=ts(jnj)         ###转化为时序格式 
plot(jnj,type="o") ###画时序图
plot(log(jnj),type="o")
xt=log(jnj)
par(mfrow=c(2,2))
acf(xt,main="x")
dxt=diff(xt)
acf(dxt,main="dx")

dsxt=diff.ts(xt,lag=4)
acf(dsxt)
dxds=diff.ts(dxt,lag=4)
acf(dxds)

#######Example 2.3
(airmod=arima(xt,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4),include.mean=F,method="ML"))
Box.test(airmod$resid,12,type="Ljung")
1-pchisq(10.748,10)

(airmod2=arima(xt[1:76],order=c(0,1,1),seasonal=list(order=c(0,1,1),period=4)))
airpred=predict(airmod2,n.ahead=8)

U=airpred$pred+1.916*airpred$se
L=airpred$pred-1.96*airpred$se
par(mfrow=c(1,1))
ts.plot(xt,airpred$pred,col=1:2,xlim=c(70,84),ylim=c(1.5,3))
lines(U,col="blue",lty="dashed")
lines(L,col="blue",lty="dashed")

###Example 2.4
decile=read.table(file="m-decile1510.txt")
Rt=ts(decile[,2])
par(mfrow=c(2,2))
plot(Rt)
acf(Rt)
pacf(Rt)

(fit1=arima(Rt,order=c(1,0,0),seasonal=list(order=c(1,0,1),period=12)))
(n=length(Rt))
(m=n/12)
jan=rep(c(1,0,0,0,0,0,0,0,0,0,0,0),m)

fit3=arima(Rt,xreg=jan,order=c(0,0,1))
e3=fit3$resid
ts.plot(e3)
acf(e3)

############################带时序误差的回归模型
dat=read.table(file="w-gs1n36299.txt")
r1=ts(dat[,1])
r3=ts(dat[,2])
plot(r1,col="red",type="l")
lines(r3,col="blue",lty="dashed")
plot(r1,r3)
lm.r=lm(r3~r1) ###直接做回归
summary(lm.r)
plot(lm.r$resid,type="l")
acf(lm.r$resid) 

c1=diff(r1)
c3=diff(r3)

par(mfrow=c(1,2))
plot(r1,r3)
plot(c1,c3)

lm.c=lm(c3~c1)
summary(lm.c)
plot(lm.c$resid,type="l")
acf(lm.c$resid)

fit3=arima(c3,xreg=c1,order=c(0,0,1))
acf(fit3$resid)
#################################长记忆模型

